Applications of SOS Programming in Flowpipe Construction

Flowpipe construction consists of under or over-approximating the sets of states reachable by dynamical systems. Recently a method has been developed for the class of polynomial ODEs with uncertain initial states (see [1], abbreviated XFZ18under). This method consists of reducing the Hamilton-Jacobi-Bellman equation to a hierarchy of semidefinite programs.

In this notebook we consider the problem of approximating the flowpipe of a system of polynomial ODEs using XFZ18under. This is a Julia implementation that relies on the JuMP ecosystem (JuMP, PolyJuMP, SumOfSquares, MathProgInterface, MathOptInterfaceMosek) and the JuliaAlgebra ecosystem (MultivariatePolynomials, DynamicPolynomials). The implementation is evaluated on a set of standard benchmarks from formal verification and control engineering domains.


References:

Quick links to documentation:

Van-der-Pol system

Consider the following van-der-Pol system:

$$ \dot{x}_1 = x_2 \\ \dot{x}_2 = -0.2x_1 + x_2 - 0.2x_1^2 x_2 $$

In [1]:
using MultivariatePolynomials,
      JuMP,
      PolyJuMP,
      SumOfSquares,
      DynamicPolynomials,
      MathOptInterfaceMosek,
      MathematicalSystems

const  = differentiate
Out[1]:
differentiate (generic function with 18 methods)
In [2]:
# symbolic variables
@polyvar x₁ x₂ t

# time duration (scaled, see dynamics below)
T = 1.0 

# dynamics
f = 2 * [x₂, -0.2*x₁ + x₂ - 0.2*x₁^2*x₂] 

# set of initial states X₀ = {x: V₀(x) <= 0}
V₀ = x₁^2 + x₂^2 - 0.25

# constraints Y = {x: g(x) >= 0} compact search space Y x [0, T]
g = 9 - x₁^2 - x₂^2

# degree of the relaxation
k = 4

# monomial vector up to order k, 0 <= sum_i alpha_i <= k, if alpha_i is the exponent of x_i
X = monomials([x₁, x₂], 0:k)
XT = monomials([x₁, x₂, t], 0:k)

# create a SOS JuMP model to solve with Mosek
model = SOSModel(with_optimizer(MosekOptimizer))

# add unknown Φ to the model
@variable(model, Φ, Poly(XT))

# jacobian
∂t = α -> (α, t)
∂xf = α -> (α, x₁) * f[1] + (α, x₂) * f[2] 
 = ∂t(Φ) + ∂xf(Φ)

# Φ(x, t) at time 0
Φ₀ = subs(Φ, t => 0.)

# scalar variable
@variable(model, ϵ)

@variable(model, s₁, SOSPoly(XT))
@variable(model, s₂, SOSPoly(XT))
@variable(model, s₄, SOSPoly(XT))
@variable(model, s₅, SOSPoly(XT))

@variable(model, s₇, SOSPoly(X))
@variable(model, s₉, SOSPoly(X))

@constraint(model, ϵ >= 0.)
@constraint(model,  - s₁*t*(T-t) - s₂*g  SOSCone())
@constraint(model, ϵ -  - s₄*t*(T-t) - s₅*g  SOSCone())
@constraint(model, Φ₀ - V₀ - s₇*g  SOSCone())
@constraint(model, ϵ + V₀ - Φ₀ - s₉*g  SOSCone())

@objective(model, Min, ϵ)
Out[2]:
$$ ϵ $$
In [3]:
optimize!(model)
MOSEK warning 705: #1 (nearly) zero elements are specified in sparse row ''(7043) of matrix 'A'.
MOSEK warning 705: #1 (nearly) zero elements are specified in sparse row ''(7044) of matrix 'A'.
MOSEK warning 705: #1 (nearly) zero elements are specified in sparse row ''(7045) of matrix 'A'.
MOSEK warning 705: #1 (nearly) zero elements are specified in sparse row ''(7046) of matrix 'A'.
MOSEK warning 705: #2 (nearly) zero elements are specified in sparse row ''(7047) of matrix 'A'.
MOSEK warning 705: #2 (nearly) zero elements are specified in sparse row ''(7048) of matrix 'A'.
MOSEK warning 705: #2 (nearly) zero elements are specified in sparse row ''(7049) of matrix 'A'.
MOSEK warning 705: #3 (nearly) zero elements are specified in sparse row ''(7050) of matrix 'A'.
MOSEK warning 705: #3 (nearly) zero elements are specified in sparse row ''(7051) of matrix 'A'.
MOSEK warning 705: #4 (nearly) zero elements are specified in sparse row ''(7052) of matrix 'A'.
Warning number 705 is disabled.
Problem
  Name                   :                 
  Objective sense        : min             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 7119            
  Cones                  : 0               
  Scalar variables       : 6450            
  Matrix variables       : 10              
  Integer variables      : 0               

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 390
Eliminator terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries                  : 2                 time                   : 0.00            
Lin. dep.  - tries                  : 1                 time                   : 0.00            
Lin. dep.  - number                 : 0               
Presolve terminated. Time: 0.02    
GP based matrix reordering started.
GP based matrix reordering terminated.
Problem
  Name                   :                 
  Objective sense        : min             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 7119            
  Cones                  : 0               
  Scalar variables       : 6450            
  Matrix variables       : 10              
  Integer variables      : 0               

Optimizer  - threads                : 8               
Optimizer  - solved problem         : the primal      
Optimizer  - Constraints            : 6728
Optimizer  - Cones                  : 1
Optimizer  - Scalar variables       : 6061              conic                  : 6060            
Optimizer  - Semi-definite variables: 10                scalarized             : 6414            
Factor     - setup time             : 2.19              dense det. time        : 0.31            
Factor     - ML order time          : 0.22              GP order time          : 1.18            
Factor     - nonzeros before factor : 3.43e+06          after factor           : 4.44e+06        
Factor     - dense dim.             : 2                 flops                  : 4.35e+09        
Factor     - GP saved nzs           : 2.91e+06          GP saved flops         : 5.50e+09        
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  
0   1.0e+00  4.0e+00  1.0e+00  0.00e+00   0.000000000e+00   0.000000000e+00   1.0e+00  2.22  
1   5.7e-01  2.3e+00  8.7e-01  1.47e+00   6.331642200e-01   1.431120861e-02   5.7e-01  2.60  
2   1.1e-01  4.6e-01  9.0e-01  1.78e+00   5.062461633e-01   3.303437431e-02   1.1e-01  2.98  
3   2.3e-02  9.1e-02  6.6e-01  2.03e+00   2.292668064e-01   1.650739108e-01   2.3e-02  3.40  
4   4.9e-03  1.9e-02  2.4e-01  1.04e+00   5.760389317e-01   5.629009543e-01   4.9e-03  3.76  
5   1.5e-03  6.1e-03  8.4e-02  3.85e-01   8.411574667e-01   8.363075016e-01   1.5e-03  3.99  
6   4.1e-04  1.6e-03  2.4e-02  1.79e-01   1.154237607e+00   1.153934851e+00   4.1e-04  4.21  
7   1.3e-04  5.0e-04  7.8e-03  1.23e-01   1.460241945e+00   1.461710063e+00   1.3e-04  4.42  
8   3.3e-05  1.3e-04  2.0e-03  3.20e-03   1.879243467e+00   1.882100664e+00   3.3e-05  4.67  
9   1.0e-05  4.2e-05  6.4e-04  1.99e-02   2.316545761e+00   2.319942755e+00   1.0e-05  4.89  
10  3.1e-06  1.2e-05  1.9e-04  9.89e-02   2.881236163e+00   2.884922049e+00   3.1e-06  5.11  
11  1.1e-06  4.2e-06  7.1e-05  6.73e-02   3.551447326e+00   3.554691958e+00   1.0e-06  5.32  
12  4.7e-07  1.8e-06  3.6e-05  3.06e-01   4.012134792e+00   4.014602498e+00   4.6e-07  5.54  
13  1.8e-07  4.8e-07  1.7e-05  5.91e-01   4.541036416e+00   4.541835263e+00   1.2e-07  5.76  
14  1.4e-07  6.7e-08  5.3e-06  8.70e-01   4.762159761e+00   4.762315429e+00   1.8e-08  5.98  
15  5.4e-08  1.6e-08  2.5e-06  9.82e-01   4.796421007e+00   4.796458861e+00   4.1e-09  6.66  
16  4.4e-08  4.4e-09  1.3e-06  9.95e-01   4.804848574e+00   4.804859200e+00   1.1e-09  6.91  
17  4.4e-08  4.4e-09  1.3e-06  9.98e-01   4.804848574e+00   4.804859200e+00   1.1e-09  7.17  
Optimizer terminated. Time: 7.26    

In [4]:
println("Relaxation order : k = $k")
println("JuMP.termination_status(model) = ", JuMP.termination_status(model))
println("JuMP.primal_status(model) = ", JuMP.primal_status(model))
println("JuMP.dual_status(model) = ", JuMP.dual_status(model))
println("JuMP.objective_bound(model) = ", JuMP.objective_bound(model))
println("JuMP.objective_value(model) = ", JuMP.objective_value(model))
Relaxation order : k = 4
JuMP.termination_status(model) = SlowProgress
JuMP.primal_status(model) = FeasiblePoint
JuMP.dual_status(model) = FeasiblePoint
JuMP.objective_bound(model) = 0.0
JuMP.objective_value(model) = 4.804848573577511
In [5]:
# Recovering the solution:
ϵopt = JuMP.objective_value(model)

# Punder <= 0
Punder = subs(JuMP.value(model[:Φ]), t => T)
Out[5]:
-0.03914340589015114x₁^{4} + 0.0589251230593361x₁^{3}x₂ - 0.025813099537286843x₁^{2}x₂^{2} + 0.004310739496267609x₁x₂^{3} - 0.013301479232761591x₂^{4} - 5.132051053473244e-19x₁^{3} - 4.3746579200085004e-18x₁^{2}x₂ + 6.570198634009879e-18x₁x₂^{2} - 5.105246868450041e-18x₂^{3} + 0.7251205878772357x₁^{2} - 1.1289850430578001x₁x₂ + 0.7193493480163449x₂^{2} - 3.456597206203931e-17x₁ + 5.926177065379741e-17x₂ + 5.947011573240102
In [6]:
# Pover <= 0
Pover = subs(JuMP.value(model[:Φ]), t => T) - ϵopt * (T+1)
Out[6]:
-0.03914340589015114x₁^{4} + 0.0589251230593361x₁^{3}x₂ - 0.025813099537286843x₁^{2}x₂^{2} + 0.004310739496267609x₁x₂^{3} - 0.013301479232761591x₂^{4} - 5.132051053473244e-19x₁^{3} - 4.3746579200085004e-18x₁^{2}x₂ + 6.570198634009879e-18x₁x₂^{2} - 5.105246868450041e-18x₂^{3} + 0.7251205878772357x₁^{2} - 1.1289850430578001x₁x₂ + 0.7193493480163449x₂^{2} - 3.456597206203931e-17x₁ + 5.926177065379741e-17x₂ - 3.662685573914919
In [7]:
using ImplicitEquations, Plots
gr() # or pyplot()

G = plot()

_Punder(x, y) = sum([Punder.a[i]*x^exponents(mi)[1]* y^exponents(mi)[2] for (i, mi) in enumerate(Punder.x)])
plot!(G, _Punder  0., xlims=(-10, 10), ylims=(-10, 10), color="red")

_Pover(x, y) = sum([Pover.a[i]*x^exponents(mi)[1] * y^exponents(mi)[2] for (i, mi) in enumerate(Pover.x)])
plot!(G, _Pover  0., xlims=(-10, 10), ylims=(-10, 10), color="blue")

G
Out[7]:
-10 -5 0 5 10 -10 -5 0 5 10
In [8]:
G = plot()

_Punder(x, y) = sum([Punder.a[i]*x^exponents(mi)[1]* y^exponents(mi)[2] for (i, mi) in enumerate(Punder.x)])
Gu = plot(_Punder  0., xlims=(-8, 8), ylims=(-8, 8))

_Pover(x, y) = sum([Pover.a[i]*x^exponents(mi)[1] * y^exponents(mi)[2] for (i, mi) in enumerate(Pover.x)])
Go = plot(_Pover  0, xlims=(-8, 8), ylims=(-8, 8))

plot(Gu, Go)
Out[8]:
-6 -3 0 3 6 -6 -3 0 3 6 -6 -3 0 3 6 -6 -3 0 3 6

Comparison with a MATLAB/YALMIP implementation

Package k Constraints Scalar variables Matrix variables Time(s)
JuMP 2 1017 803 10 < 1
YALMIP 2 152 63 10 < 1
JuMP 3 2871 2471 10 ~ 1
YALMIP 3 254 121 10 ~ 1
JuMP 4 7119 6450 10 6.22
YALMIP 4 394 206 10 1.18